Multifractality can be a universal signature of phase transitions 
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Macroscopic systems often display phase transitions where certain physical quantities are singular or self- 
similar at different (spatial) scales. Such properties of systems are currently characterized by some order 
parameters and a few critical exponents. Nevertheless, recent studies show that the multifractality, where a 
large number of exponents are needed to quantify systems, appears in many complex systems displaying self- 
similarity. Here we propose a general approach and show that the appearance of the multifractality of an order 
parameter related quantity is the signature of a physical system transiting from one phase to another. The distri- 
bution of this quantity obtained within suitable (time) scales satisfies a g-Gaussian distribution plus a possible 
Cauchy distributed background. At the critical point the q-Gaussian shifts between Gaussian type with nar- 
row tails and Levy type with fat tails. Our results suggest that the Tsallis g-statistics, besides the conventional 
Boltzmann-Gibbs statistics, may play an important role during phase transitions. 
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Phase transitions are ubiquitous in physical systems. A first 
general theory of phase transitions is the mean field theory, 
from which in the vicinity of the critical point the free en- 
ergy can be expanded in a power series of a universal quantity 
called "order parameter" yj. This quantity fluctuates around 
zero in a "disordered" phase above the critical point, and is 
non-zero in an "ordered" phase below the critical point rep- 
resenting the broken symmetry. Unfortunately the predictions 
from this theory are often not correct when comparing with 
experiments. To remedy this, much progress has been made in 
the theory of phase transitions during last fifty years with the 
landmark of the discovery of some universal scaling laws and 
the introduction of the renormalization group method |l2)-Q| . 
With them the seemingly complex phase transitions in a va- 
riety of systems can be quantified by a few critical exponents 
which represent the self-similar characteristics of systems dur- 
ing transitions. Recently, such ideas of statistical physics have 
been applied to many complex systems in various fields and 
have extensively expanded the understanding to these sys- 
tems IJSll- Yet, these studies indicate that in many cases a few 
scaling exponents representing the self-similarity fail to quan- 
tify certain property of systems. Instead, one needs a large 
number (spectrum) of scaling exponents to clarify the char- 
acteristic of systems which represents a higher level of com- 
plexity. Examples include the processes of diffusion-limited 
aggregation (DLA) M, turbulence and chaos llllitl. Human 
heartbeat dynamics Q, climate change jlOll and many finan- 
cial quantities such as stock prices ^. This characteristic 
has been called "multifractality". A natural question is thus 
whether such behaviour can be seen in phase transitions. Dur- 
ing last two decades people have verified the multifractality of 
certain quantities on some specific transitions such as random 
field transitions llllll and Anderson transitions at the critical 
point 1I12I - U5I1 . However, a general interpretation is still ab- 
sent on how the multifractality appears with phase transitions. 
Especially, is the multifractality a universal property of phase 
transitions as those scaling laws? To address these questions, 
we argue that one should investigate the property of systems 
which appears universally in phase transitions. The only ideal 
quantity to our knowledge is the order parameter. 

Another important progress recently in statistical physics 



eralization of the conventional Boltzmann-Gibbs statistics. 
Specifically, when the effective number of degrees of freedom 
is small, one obtains the Tsallis (/-statistics following the same 
arguments to obtain the Boltzmann-Gibbs statistics lll8l[l9ll . 
In this frame some variables such as the entropy and the en- 
ergy, may be non-extensive when systems have long-range 
correlations pol]. It is well-known that the correlation length 
is infinite at the critical point of a continuous transition. Thus, 
systems near the critical point are ideal candidates to observe 
such non-extensive behaviour Just as the limit distribution of 
the sum of independent and identically distributed variables is 
Gaussian or Levy-stable distribution, in g-statistics where the 
long-range correlation is present we observe a q-Gaussian dis- 
tribution 1I21I1 . The distributions related to Tsallis q-statistics 
have been observed in many experiments, such as cold atoms 
in dissipative optical lattices 1I22I1 . superdiffusion Ii23il , solar 
plasma dvnami cs |l24ll . spin glass relaxation luSll . tissue radi- 
ation response ll26ll . and financial signals IdI uTIi . 

Here we show that, at temperatures near the critical point 
the distribution of an order parameter related quantity satisfies 
a q-Gaussian distribution plus a possible Cauchy background. 
At the critical point, the distribution shifts between the Levy 
regime and the Gaussian regime 11281 [2911 . triggering the 
multifractal behaviour and signalling the phase transition. 

Distribution and tlie measuring quantity 

The q-Gaussian distribution of interest is a symmetric distri- 
bution and has the following form: 



fg{D) = Fo • [1 + (9 - 1){D - Z?o)V2^' 
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where tr is a scale parameter related to the variance of the 
distribution. Do indicates the peak position, and Pq is a nor- 
malization parameter Letting q ^ 1 one recovers the Gaus- 
sian distribution. When q e [1,5/3), the distribution is in 
the Gaussian regime. The signal is monofractal and can be 
described by a single exponent. The distribution falls into 
the Levy regime with infinite variance when q e (5/3, 3), 
and the critical value is qc = 5/3 111711 . In this regime 
fq(D) decays asymptotically as that of a Levy-stable distri- 
bution. Specifically the Cauchy distribution coiTesponds to 
f2{D) = j/{n[{D - D^f + 7^]} where 7 = V^a is a 
scale parameter It has been shown that the signal with the 



q-Gaussian distribution in the Levy regime is multifractal 11281 - 
[3011 . The multifractal behaviour is manifested by two attrac- 
tors at (a, /(a)) == (0,0) and ((g - l)/(3 - q),l), where 
a is the singularity strength and f{a) is the singularity spec- 
trum 12911 . Further, the multifractality may survive even when 
the time correlation is destroyed by shuffling the data 1291. 

We hypothesize a universal measuring quantity which pos- 
sibly satisfies a q-Gaussian distribution as follows. It con- 
tains only the order parameter and does not depend on details 
of a specific system, thus should be a dimensionless quan- 
tity. We further utilize the "universal" self-similar property 
of certain physical quantities of the system near the phase 
transition. This implies that such property does not change 
when the measuring scale alters. Combining these ideas we 
construct the desired quantity as the ratio of two order pa- 
rameters measured at different spatial scales. Its distribution 
P{D) = P{mi/mi,) where mi and mi, are order parameters 
m measured at scale 1 and scale b > 1, respectively. 

When far away from the critical point the distribution of an 
order parameter should be Gaussian at different spatial scales. 
Its mean is zero in a disordered phase and non-zero in an or- 
dered phase. Thus in the disordered phase P{D) satisfies a 
Cauchy distribution. In the ordered phase, we mark the mean 
and the fluctuations of the order parameter at the scale "1" 
(scale "fe") as xq and Sx{t) [j/o and Sy{t)], respectively. For 
sufficiently large systems xq ^ Sx{t) and yo ^ Sy{t). Thus 
the ratio D is: 



Xq + Sx{t) 
yo + 6y{t) 



Xq^ y0^x{t) - XoSy{t) 
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xoSyitf - yQSxit)Sy{t) 
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From equation dU we find that in the ordered phase P{D) to 
the first order of (5 is a Gaussian. Therefore, in two extreme 
cases, P{D) satisfies a g-Gaussian distribution. We next 
investigate how P{D) behaves when the system is near the 
critical point from some tunable examples. 

Example tunable models 

We now focus on some tunable examples of finite size spin 
systems to verify our hypothesis. For all models we obtain 
the data through Monte Carlo simulations and apply peri- 
odic boundary conditions to systems of different sizes. When 
the size L — >^ 00 one obtains the property of a macroscopic 
system. The system at the scale "1" is the original spin 
system, and the block spin system at the scale "6" is con- 
structed using the standard coarse-graining operation in sta- 
tistical physics [2-4] (see the Appendix[A]i. 

We investigate two well-known models: (1) The 2- 
dimensional nearest-neighbour Ising model with the Hamil- 
tonian H ~ ^JJ2(ij)^i^j^ where (ij) represents nearest 
neighbours, ferromagnetic parameter J = 1 and the spin 
Si = ±1. This model has a second order phase transition 
at Tc = 2 J/[ln(l + \/2)]- The temperature is in units of 1/k^ 
where fee is the Boltzmann constant. Its order parameter is the 
magnetization (per site) and we obtain the time series using 



flie Wolff algorithm Jilt]. The ratio D{t) = m(l, t)/m(6, t) 
where m{b, t) and m(l, t) are order parameters measured at 
the same Monte Carlo step (MCS) "t". (2) The 3-dimensional 
nearest-neighbour Ising glass model which has a second or- 
der phase transition with very slow dynamics near the crit- 
ical point iH mt]. Its Tc ^ 0.95 and the Hamiltonian 
H = —J2(ii) JijSiSj where we set Jij satisfying a Gaus- 
sian distribution with zero mean and unit variance, and the 
spin Si = ±1. Each configuration of Jy is one sample. At 
each temperature after the equilibrium we simulate 150 sam- 
ples with at least two million MCS per sample. The order 

parameter is the spin overlap rn = I ^^ s^ ' s^ ') /N where 

{Si } is a replica of {s^ }, N is the number of spins. The 
order parameter time series are obtained using the Metropolis 
algorithm iQ. 

We also provide another model as a representing scenario 
for general cases of phase transitions. The example is the 2- 
dimensional nearest-neighbour n-state Potts model where n is 
a positive integer and each spin takes values 0,l,...,n — 1. 
This model has a second order transition for n < 4 and a first 
order transition for ?i > 4 Il35ll . Specifically, when n = 2 
it returns to the Ising model. Its T^ = J/[ln(l + \/n)], the 
Hamiltonian H = ^J^Uj) Ss,.sj and we set J = 1. The 
order parameter is nii^a = {n x {6{si,a))T — ^)/{n— 1) 
where i is the index of a spin, a is the value that the spin 
can take and (. . .)t indicates a thermal average Ii36ll . The 
order parameter series are obtained using the Metropolis al- 
gorithm 1341] . Due to the spin symmetry, we obtain P{D) uti- 
lizing all Dij{t) — mi^Q,(l, t)/'mj,a{b, t) where j is the index 
of the block spin, i is the index of any original spin within the 
block spin j. We have fixed the value of a in the calculation. 

In the simulations we set the rescaling factor 6 = 4 for 
the Ising model and the Potts model, and 6 = 2 for the glass 
model. 

P{D) near the critical point 

Results of P{D) for the Ising model near the critical point are 
shown in Fig.[T^ and in the Appendix [B] When approaching 
the critical point from the ordered phase (T < T^), a Gaus- 
sian fit to P{D) gradually does not work well, and it is worse 
when T > Tc. In these situations the testing Gaussian fit un- 
derestimates P{D) at the peak, overestimates P{D) on two 
shoulders, and diminishes faster than the exponential decay 
compared to the decay of P{D) which contains power-law 
tails. Instead, a g-Gaussian fit coincides very well with the 
central part of P{D) at all temperatures. We observe an in- 
creasing value of the parameter q with increasing temperature. 
At q = 5/3 the q-Gaussian enters Levy regime with infinite 
variance, and we take this point as our critical temperature 
T^„ for the finite system with the size L. Thus from this spe- 
cific example, we find that near the critical point the order 
parameter ratio D is our desired quantity. 

When the ratio D is very far away from the peak Dq, we 
find that P{D) separates from such q-Gaussian fit. At all tem- 
peratures we observe that fat tails of P{D) decay as D~^, 
suggesting possible "universal" origin. Since in the disor- 
dered phase P{D) follows a Cauchy distribution, we fit these 
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FIG. 1: Distributions of tlie order parameter ratio P{D). a, At dif- 
ferent temperatures for the 2-dimensional Ising model with the size 
L = 124. b, At the critical temperature T^g for the 3-dimensional 
Ising glass model with the size L = 12, and the 2-dimensional 5- 
state Potts model with the size L = 48 and the thermal averaging 
interval A'tw ~ 80, 000. The lengths of the order parameter ratio 
for three models are 8-16 million, 150 samples each with 2 million 
points, and 5.8 million, respectively. At all temperatures the central 
part of P{D) can be well fitted by a g-Gaussian distribution, and the 
fat tails can be well fitted by a Cauchy distribution. Summation of 
the two provides a good fit at all ranges and at all temperatures. A 
Gaussian fit (dashed line) to the central part does not work when T 
approaches T^^ and when T > T^q. More details of the tails for the 
same distributions see the Appendix O 



fat tails with this distribution and take them as background. 
The Cauchy fits work well even when P{D) ~ 10^^ (see 
the Appendix FBIi. As shown in Fig. [T^, the summation of a 
g-Gaussian distribution and a suitable Cauchy background at 
each temperature of interest provides good fit to P{D) in all 
ranges of D, i.e., P{D) = {I - p) /,(£>) + p f2{D) where 
/2(D) is due to the Cauchy tails and p is the probability of 
the Cauchy contribution. This implies that the central part 
and tails part of P{D) may be independent. By reversing our 
procedure we can deduce the origin of the Cauchy part (see 
the Appendix fell. We find that the distribution of m(l,t) or 
m{h, t) which contributes to the Cauchy part outside the cross 
points with the central part of P{D) achieves a local maxi- 
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FIG. 2: Characteristics of P{D) for the 5-state Potts model measured 
with different thermal averaging interval 7Vt„. a. Asymmetric and 
symmetric distributions with different A'^t„ for the 5-state Potts model 
at T = 0.8495 (g = 1.3). The system size is L = 48 and the red 
lines are g-Gaussian fits, b. Dependence of different parameters on 
the value of A'tw where we fix the temperature T = 0.8515. 



mum and is symmetric about zero, indicating a disorder-like 
behaviour It decays with approximately exponential tails for 
large values of m and decays faster for the larger system. 

Similar considerations can be done on both the Ising glass 
model and the Potts model. We find that the ratio distributions 
of them both share the similar behaviour with that of the Ising 
model. The Gaussian distribution could not fit the central part 
of P{D) at temperatures near the critical temperature T^ . 
This is true even for the 5-state Potts model which has a first 
order phase transition. In Fig. [TJ) we show the fits at T^ for 
both models. When approaching T^ we observe Cauchy 
distributed fat tails where the ratio D is far away from the 
peak Do- Further, combining a g-Gaussian distribution and 
a suitable Cauchy background provides a good fit to P{D) 
at all ranges of D. Nevertheless, for the Potts model, one 
has to choose a suitable thermal averaging interval iVt„ (in 
units of MCS) for the order parameter For the 5-state Potts 
model with the size L = 48 we take TVtw = 80, 000 at all 
temperatures. 

General scenario 

As shown in Fig.|2^, P{D) taken with randomly chosen A^tw, 
e.g., iVtw = 2, 000 for the 5-state Potts system with the size 
L ~ 48, may be asymmetric thus would not follow a q- 
Gaussian distribution. With the Potts model as an example, 
here we propose an approach which can be applied to general 
physical systems. 

First, from the construction of our measuring quantity we 
hypothesize that the possible "universal" g-Gaussian distribu- 
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FIG. 3: Characteristics of tlie Caucliy noise strength in P{D) within 
suitable range of A'tw The Cauchy noise strength is proportional 
to the ap of the noise where D is far away from the peak Do and is 
invariant within suitable range of A'tw Here we show the dependence 
of the ap of the Cauchy noise on the system size L for fixed q, and its 
dependence on the value of q for fixed L = 124 and Nt„ — 15, 000. 



tion is relevant to the self-similarity of certain physical quanti- 
ties near the transition. To verify this, we investigate the self- 
similarity of the distribution P(m) of the order parameter m 
at different spatial scales and find that it is broken for small 
values of A'tw (see the Appendix[A]i. Correspondingly we ob- 
tain asymmetric P{D). The self-similarity of P{m) becomes 
better for larger values of A'tw. When A'tw = 150, 000 for 
the 5-state Potts model with the size L = 96 the distributions 
P{m) at two scales are almost identical after the rescaling. 
Correspondingly P{D) becomes symmetric. 

We thus next quantify how P{D) responds when the value 
of A'tw alters. To do this, we calculate the skewness of the cen- 
tral part of P{D). As examples shown in Fig.|2^, the consid- 
ered region is within two dashed lines. We fix the temperature 
and study how the skewness changes with different values of 
A'tw (see the Appendix FDI for relevant details). For different 
sizes of systems we find that the skewness is closer to zero 
for larger value of iVtw. For sufficiently large A'tw > ^tw" the 
skewness stays at around zero, while P{D) is symmetric and 
its central part follows a g-Gaussian distribution. The value of 
A'tw'" seems larger for larger size of systems and it is compara- 
ble with the characteristic time scale of the order parameter m 
(see the AppendixFEJi. Specifically, for the 2-state Potts model 
A',™" is around the number of spins (see the Appendix FDT i. 
thus the results of the 2-state Potts model is consistent with 
that of the Ising model. Starting from A't™", some specific 
short-range properties of the Potts model are lost and the q 
value from a q-Gaussian fit to P{D) keeps constant in a broad 
range of A'tw (see Fig. ^p). Thus we consider the value of q 
obtained with A'tw chosen in such range as the suitable value 
of q. When A'tw further increases, however the value of q will 
increase. In this situation the order parameter itself is gradu- 
ally not a good quantity to quantify the system. The above are 
remarks how we should choose the order parameters. 
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FIG. 4: Critical behaviour for all three models. The systems enter 
the Levy regime when q > qc = 5/3. When q < 5/3, the systems 
are in the Gaussian regime. For better view we have made some scale 
transfomiations: t -)> 5i for the 5-state Potts (with A'tw = 80, 000) 
and t — i- t/5 for the glass. Inset: Finite size scaling for all three 
models. The dashed lines are fitting lines. We show how the size- 
dependent critical temperature T^ 



We mark t = (T^ - %,q)/T^. 



depends on the system size L. 
'here Tc,q = T^. We have made 
some scale transformations: f — >■ 4t and L""' — >• AL~^'' for the 
5-state Potts model; t ->• i/20 and L^"" -s> L^"" /IQ for the Ising 
glass model. We observe power-law dependence for all models. 



Finally we evaluate the effect of the Cauchy noise which 
appears for all models. It is pronounced where D is far away 
from Do and its contribution is p/2 (I?) ~ ap/{D — Daf'. We 
thus consider ap as a reliable quantity to measure the strength 
of the Cauchy noise. In the suitable range of A'tw, the crp of the 
Cauchy noise is invariant (see Fig.|2j5 and the Appendix [D). 
We then study how it depends on the system size L for fixed 
q and suitable A'tw. As examples shown in Fig. [3] for both the 
2-state Potts model from our general approach and the Ising 
model we find that it decays in power-law with the system 
size. This implies that for a macroscopic system the Cauchy 
contribution may go to zero and P{D) is (/-Gaussian in all 
ranges of D. We also study how the ap of the Cauchy noise 
depends on the parameter q for fixed L and A'tw. By com- 
paring the results of two equivalent models mentioned above, 
we find that P{D) from the former contains weaker Cauchy 
noise. Further, we also obtain better g-Gaussian fit to P{D) 
of the Potts model which is constructed with fewer points (see 
Fig.ffl). 

Critical behaviour and multifractality 

We now subtract the Cauchy noise and investigate the rest of 
P{D), whose behaviour is related to macroscopic systems. 
In Fig. |4] we show for all three models how the value of q 
from a q-Gaussian fit depends on the reduced temperature 
t ^ [T — Tc)/Tc, where T^ is obtained from the conventional 
methods. For each finite system with the size L, we search for 
its corresponding critical temperature T^ . We employ this on 
different sizes of systems and all models. From the finite size 
scaling we have T^ = aL~^'' + T^.q. We take the values of 
Tc ,j as those from conventional methods and find good power- 




" a = 0.00987 



■a = 0.0109 




10' 10' 10^ 10' 10' 



^0.5 



•-• 5 = 1 .97 
■-■ 5 = 1 .92 
♦-♦ 5 = 1 .83 
*r^5' = 1-75 


^^ 0.8 

0.6 

<0.4 

0.2 



1 


?« 


^ 


fe**= 


"^ 












■ 


2 


1.4 


1.6 


1.8 ; 





0.2 



0.4 



0.6 



a 



FIG. 5: Multifractal behaviour of the Ising model near the critical 
point. The system size is i = 124 and we have applied the MFDFA. 
The data have been shuffled before the calculation and the length of 
the order parameter ratio is 2 million, a, The fluctuation function 
Fr(n) versus observing window n at different temperatures. Both 
the parameters q and a are obtained from the g-Gaussian fit to the 
central part of P{D). Fr{n) are calculated in the ranges of r G 
[-4, 4] and n G (20, 327, 680). b, The singularity spectrum f{a) 
versus the singularity strength a for different values of q. The fitting 
range of n to obtain a and f{a) is (20, 10,240). Inset: the range of 
the spectrum Aa — Omax — amin for different values of q. 



law behaviours between T^ — T^^g and L (see inset of Fig. |4]i. 
The critical exponents i^q for the Ising model, the Ising glass 
model and the 5-state Potts model are 1.16, 1.22 and 1.61, re- 
spectively. Thus the critical temperatures from our method is 
consistent with those from conventional considerations. 

We next quantify the multifractality in signals u sing the 
multifractal detrended fluctuation analysis (MFDFA) ll30ll37l - 
I39II (see the Appendix [Fji, which is considered to be a re- 
liable method to quantify the multifractal characteristic in a 
nonstationary series ^. With this method we can calculate 
the r-th order fluctuation function Fr{n) ^ n^^'^\ where n 
is the size of the observing window. Varying h{r) indicates 
the multifractality. We also calculate widely-used the singu- 
larity strength a and the singularity spectrum /(a). When 
{a, f{a)) spreads from one point to a variety of points in the 
X-Y plane, the multifractal behaviour appears. 

In Fig. |5] with the Ising model as an example we show the 
multifractal behaviour of the shuffled series of the order pa- 
rameter ratio where the time correlation is absent. Similar 
behaviour is also seen on other models (see the Appendix [Gil. 
We find that the order parameter ratio is monofractal when 



the corresponding q is far below qc = 5/3. When we ap- 
proach T^ from the ordered phase, the slope of Fr (n) gradu- 
ally shifts from small scales of n for certain values of r. Yet at 
large scales all fluctuation functions still maintain the identi- 
cal slope. At q ^ qc the multifractality appears at all scales we 
observe. Nevertheless, the range of n where the multifractal- 
ity presents shrinks when continuing increasing the value of q. 
We further show the singularity spectrum f{a) v.s. a. We find 
that the spreading range of the spectrum Aq, = a,^^^ — Omin is 
small for small q < 5/3, indicating a Gaussian-like monofrac- 
tal behaviour. However when q approaches 5/3, A^ rapidly 
increases and then slowly expands with increasing q. 

The multifractal behaviour shown above should not be 
system specific since we have considered the shuffled se- 
ries. To further verify this, we generate some g-Gaussian 
distributed artificial signals and set the same parameters 
for the artificial signals and the testing model (see the 
Appendix [Hb. We find that the results are almost identical for 
the artificial signals and the model. Thus such multifractal 
behaviour is determined by the two parameters q and <t of 
the g-Gaussian distribution. To observe the multifractality, 
it requires q > 5/3. When q > 5/3, the range where the 
multifractality presents is controlled by a. For fixed value of 
q, the larger size of system, the smaller value of cr (see the 
Appendix [H]i. Correspondingly the larger system presents 
multifractal behaviour in broader range of n. 



Conclusion 

With tunable examples we have shown that physical systems 
near the phase transition present much more complex self- 
similar behaviour represented by the multifractality of the 
ratio of two order parameters measured at different spatial 
scales. At all temperatures around the critical point the 
distribution of the ratio follows a non-extensive g-Gaussian 
distribution plus a possible Cauchy background which decays 
in power-law with the system size and may disappear for a 
macroscopic system. The q-Gaussian distribution enters the 
Levy regime at the critical point, and triggers the multifrac- 
tality at all scales we observe. We have proposed a general 
approach which relates only to the broken symmetry yielding 
zero and non-zero order parameters in different phases as 
well as the self-similar characteristics of systems near the 
transition. Thus it should be applicable to other physical 
systems. The multifractality appears for the order parameter 
TO obtained within suitable range of (time) scales where 
certain short-range properties of the specific system are lost 
and m shows good spatial self-similarity. In this situation the 
long-range correlation in the order parameter prevails and the 
ratio of order parameters follows a non-extensive q-Gaussian 
distribution. Our results suggest that the Tsallis q-statistics 
may play an important role in phase transitions. 
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Appendix 



A. Self-similarity in distributions of the order parameter 

When a system is in the vicinity of the critical point, many 
physical quantities of it display self-similar (scaling) proper- 
ties. The distribution of the order parameter also shares this 
characteristic. In Fig. |ST| we show an example for the 2- 
dimensional 5-state Potts model with the size L = 96. We 
construct the associated block spin system with size L' = 
96/4 = 24 utilizing the standard coarse-graining procedure. 
For a d-dimensional spin system with L^ lattices, we can 
transform it into a block spin system with {L/b)'^ block spins. 
Each of block spins contains b"^ spins. A coarse-graining op- 
eration is reached when we substitute each of block spins by 




0.5 1 

Order parameter m 

FIG. S 1 : Self-similarity of distributions of the order parameter for the 
5-state Potts model. The system size is i = 96 and the temperature 
is T = 0.8515. For the distribution of the block (coarse-grained) 
spin system (with L' = 96/4 — 24) we have done the following 
scale transformation: x — )■ 2;/1.215, y ^ y x 1.215. With larger 
thermal averaging interval A^tw we obtain better similarity between 
the distribution of the original spin system and that of the block spin 
system. 




FIG. S2: Cauchy distributed fat tails in distributions of the order 
parameter ratio. We show such tails at different temperatures for all 
three models. The data are the same as those we show in Fig. I of 
the paper. 



a single spin with the value determined by the majority rule. 
For example, for the Ising model the block spin is H-1 if there 
are more spins up than down, and vice versa. In particular, 
when the amount of spins up and spins down is equal, we as- 
sign the block spin h-1 or -1 randomly. In this way we have 
constructed the Ising system at the scale "&". The block spin 
system for the Potts model can be constructed similarly. 

Nevertheless, for the Potts model we do not see good scal- 
ing when the thermal averaging interval A^tw for the order pa- 
rameter m is too small (see Fig.lSlll. In this situation m con- 
tains certain short-range information specifically related to the 
Potts model. It is well-known that the critical properties of a 
physical system do not depend on the short-range details, but 
on the characteristics of long-range fluctuations. Such short- 
range information is not self-similar and may diminish with 
larger A'tw where the time correlation in the order parameter 
becomes weaker (For example, see the Appendix [E]) For 
sufficiently large value of A^tw, the distribution of the order 
parameter for the original system hi {x) and that for the block 
spin system hi,{x) become similar, i.e., hi{x) = hi,{cx)/c. 
Interestingly, when such scaling relation is effective, the dis- 
tribution of the order parameter ratio P{D) becomes symmet- 
ric and follows a (jf-Gaussian distribution, as we show in the 
paper. 



B. Cauchy noise in distributions of the order parameter ratio 

In Fig.[S2]we present the Cauchy distributed fat tails in dis- 
tributions of the order parameter ratio P{D). For the Ising 
model we observe very good fit to the tails of P{D) down to 
lO^'* — 10~^, while the length of the order parameter ratio is 



8-16 million. We also note that at high temperatures with q 
close to 2, the Cauchy background may be still pronounced. 
In an example shown in Fig. [S2] for the system at T = 2.31 
and with the size L = 124, we find that a single g-Gaussian fit 
works well for the central part, however it still underestimates 
P{D) at positions of D far away from the peak Do- 

The above behaviour has also been found in the Ising glass 
model and the Potts model. For the former at each temperature 
we have 150 samples each with 2 million points after the equi- 
librium. For the latter at each temperature the length of the or- 
der parameter ratio is 5.8 million. As examples in Fig. |S2|we 
show such behaviour at the critical point T^ which we have 
defined in the paper For both models good fits to the tails of 
P{D) down to 10^'' are observed. 



C. Origin of the Cauchy noise 

By reversing our procedure we can deduce the origin of the 
Cauchy distributed fat tails. To see this, we define a quantity 
r{D) = min[/cauchytit(-D)/P(i:>), 1], where /cauchyfit(£') is our 
Cauchy fit and P{D) is the distribution of the model from the 
Monte Carlo simulation. At each Monte Carlo step t with the 
order parameter ratio D{t) = m{l,t)/m{b,t), we take r{D) 
as the probability of this data point belonging to the Cauchy 
noise. We did this at each t and thus could obtain the cor- 
responding ?Ti(l, t) and m(6, t) which belong to the Cauchy 
noise. Such considerations can be done in all ranges of D or 
in partial range of D, e.g., we can choose the range of D which 
is outside the cross points of the central part and the Cauchy 
part of P{D). As shown in the orange and brown dashed lines 
of Fig. [S3] we find that the distribution of m{l, t) or m{b, t) 
which is associated with the Cauchy part of P{D) outside the 
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FIG. S3: Origin of the Cauchy noise as seen in the distribution P{m) 
of the order parameter m. a, For the 2-dimensional Ising model at 
T = 2.295 with the size L — 124. b, For the 2-dimensional Ising 
model at r = 2.275 with the size L = 512. c, For the 5-state 
Potts model at T = 0.851 with the size L = 48, while the thermal 
averaging interval A'^tw is 80,000. The values of q for the distributions 
of the coiTesponding order parameter ratio are 1.66, 1.66 and 1.45, 
respectively. The orange and brown dashed lines are contributed by 
the original and block spins associated with the Cauchy part of P(D) 
outside the cross points with the central part of P{D). 



cross points with the central part of P{D) achieves a local 
maximum and is symmetric about zero, then decays with ap- 
proximately exponential tails for large values of m, indicating 
a disorder-like behaviour. Further, the tails decay faster for 
the system with larger size. This implies that the fat tails in 
P{D) may go to zero when the system size L ^ cx). Such 
behaviour is also true when the distribution P{rn) of the order 
parameter is asymmetric. 



D. Asymmetric and symmetric P{D) with different Mw 

For the n-state Potts model with n > 2 we find that the 
distribution of the order parameter ratio is not symmetric 
when the thermal averaging interval A'tw is small, as shown in 
Fig. 2a of the paper. Further, such behaviour diminishes with 
larger values of A'tw To characterize this, we can calculate the 
skewness of P{D) around the peak Z?o- (In practice, the width 
of the region we choose to calculate is around llcr — 12cr, 
where a is the scale parameter of the testing g-Gaussian fit to 
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FIG. S4: Characteristics of the distributions of the order parameter 
ratio P{D) for the Potts model, a, The skewness of the central part 
of P{D) for the 5-state Potts model at temperatures with different 
values of q. We fix the system size L = 36 and the thermal averaging 
interval 7Vi„ = 16, 000. b. The dependence of the minimal Nt„ at 
which one could obtain symmetric P(D) on the Potts state number n 
at temperatures with fixed q — 1.5. We fix the system size L = 48. 
The skewness of the corresponding P{D) is within (-0.01, 0.01). 
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FIG. S5: Strength of the Cauchy noise for the 5-state Potts model. 
We fix the temperature T = 0.8515 (q ~ 1.55) and the size L = 36. 
For different values of A'^tw within the suitable range, the strength 
represented by the ap of the noise keeps constant. 



P{D).) As shown in Fig. IS4b. such asymmetry seems larger 
at lower temperature with smaller value of the parameter q. 
Further, it is also stronger for the Potts system with more pos- 
sible spin states. This is manifested in Fig. |S4b , where we 
fix the g-Gaussian fit parameter q = 1.5. We find that A^™" 
— the minimal value of the thermal averaging interval A^tw to 
obtain symmetric P{D) — increases with increasing value of 
the state number n. Specifically, we note that A"™" is con- 
sistent with the results of the Ising model which is related to 
the 2-state case of the Potts model. Since for n < 4 the sys- 
tem has a second order phase transition, such behaviour is not 
related to the order of the transition. 

For sufficiently large thermal averaging interval A'tw, P{D) 



0.5 



\ 


\^^ \ ■ 


L= 


=36 — — ^ 





2e+04 4e+04 2e+04 4e+04 6e+04 



4e+04 8e+04 

T (MCS) 



1e+05 2e+05 
T (MCS) 



FIG. S6: Average time correlation in tlie order parameter m of tlie 
5-state Potts model. We fix the temperature T — 0.8515 and vary 
the size of systems. The time delay tq corresponding to C(r) = 1/e 
is the characteristic time scale of interest. 



becomes symmetric, and it would keep symmetric when con- 
tinuing increasing the value of A'tw, as shown in Fig. 2b of the 
paper. Further, the value of the parameter q in a g-Gaussian 
fit to P{D) also keeps almost invariant in a broad range of 
Ntv,. Here we show some examples in Fig.[S5] In this suitable 
range the strength of the Cauchy noise, marked by the ap of 
the noise, also keeps constant. 



E. Time correlation In the order parameter time series 

The order parameters which we obtain through Monte 
Carlo simulations may be correlated in time. Such correla- 
tion may be stronger when the system is in the vicinity of the 
critical point. As an example, here we investigate the 5-state 
Potts model near the phase transition. We fix the temperature 
T = 0.8515 and vary the size of systems. To obtain the time 
correlation function C(r), for different sizes of systems we 
fix the thermal averaging interval A'tw = 2, 000. The order pa- 
rameter m.i^a of the Potts model depends on the spin value a 
and the spin index i. We also fix a and for each original spin i 
we calculate the time correlation Ci{T). We show the average 
C{t) of all original spins in Fig. |S6] When the time decay r 
is small C{t) decays exponentially, i.e., C{t) ^ cxp[— t/tq] 
where tq is the characteristic time scale. Comparing with the 
results we show in Fig 2b, we find that the distribution of 
the order parameter ratio P{D) becomes symmetric when the 
thermal averaging interval Nf„ is comparable or larger than 
the characteristic time scale of the order parameter m. 



F. Multifractal detrended fluctuation analysis (MFDFA) 

The multifractal detrended fluctuation analysis (MFDFA) 
is an efficient and reliable method available to quantify the 
multifractality in a non-stationary series. For a series {xi} 
with the length N, its procedure is the following; 



(1) Calculate the profile 

Y{j) 





2=1 



(x)],j = l,...,iV, 



where (...) is the mean of {xi}. 

(2) Divide Y{i) into iV„ = vi\t{N /n) non-overlapping parts 
of length 71. The maximal value of n should be smaller than 
N/A to avoid statistically unreliable results. We perform this 
step from both the beginning and the end of the signal. Thus 
we totally obtain 2iV„ parts. 

(3) For each part with index v we calculate the local trend 
yiii) in this part by a least-square fit of the series with a poly- 
nomial function, where (, is the order of the polynomial func- 
tion. We then subtract this local trend and determine the vari- 
ance: 

n 

F\v, n)^- Y,{Y[{u l)n + z] - vi{t)}\ 

n ^ — ^ 



for i^ = 1 , . . . , Nn and 

1 " 

n ^ — ^ 



i,y-N.^)n + i]^yiii)y, 



forz/ = 7V„ + l,...,2iV„. 

(4) Average over all parts we obtain the j'-th order fluctua- 
tion function: 



k I/=l 



1/r 



r/2 



Fr{n) 
For 7' = we choose 



2Ar„ 



(5) Determine the scaling of the fluctuation function: F^ (n) ^ 
jiHr) _ If i^(^j.^ jg ^ constant, the signal is monofractal; other- 
wise it is multifractal. 

The singularity spectrum can be further calculated from: 
/(a) = r[a — h{r)] + 1, where the singularity strength a = 
h{r) + rh'{r). If the signal is monofractal, we find that a = 
h{r) = const and f{a) = 1. 



G. Multifractal behaviour of the 5-state Potts model 

As another example here we investigate the multifractal be- 
haviour of the 5-state Potts model in the vicinity of the critical 
point. We set the thermal averaging interval A'tw — 80, 000 
at all temperatures. After obtaining the series of the order 
parameter ratio, we first shuffle the data and then apply the 
MFDFA. The results of the fluctuation function Fr (n) ver- 
sus the observing window n are shown in Fig. [ST] We find 
that, below the critical point (in the ordered phase) the system 
is monofractal since the slope of Fr{n) is identical for dif- 
ferent values of 7'. When approaching the critical point from 
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FIG. S7: Multifractal behaviour of the 5-state Potts model near the 
critical point. The system size is L = 48 and we have applied the 
MFDFA. The data have been shuffled before the calculation and the 
length of the order parameter ratio is 2 million. We set the thermal 
averaging interval Nt-^ — 80, 000 at all temperatures. We show the 
fluctuation function Fr{n) versus observing window n at different 
temperatures. Both the parameters q and a are obtained from the 
g-Gaussian fit to the central part of P{D). 



the region where the multifractal behaviour presents shrinks. 
Such behaviour is similar to that of the Ising model. Never- 
theless, for the Potts model we find that the scale parameter a 
of the q-Gaussian fit is increasing with increasing temperature 
and fixed A'tw, while it is not sensitive to the temperature for 
the Ising model if the measured q is far below 2. 

H. Characteristics of the multifractal behaviour 

Here we investigate the characteristics of the multifractal 
behaviour we show in the paper To do this, we generate some 
q-Gaussian distributed artificial signals and we set the same 
parameters for the artificial signals and for the testing model 
— the Ising model. For the model we have first shuffled the 
data thus eliminated the dynamics in signals. As shown in 
Fig. [S8]we find that the results are almost identical for the 
artificial signals and the model when the parameters q and 
a are the same. This verifies that the multifractal behaviour 
we present in the paper is not a system specific behaviour but 
should be applicable to other physical systems. 
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FIG. S8: Multifractal behaviour of the Ising model compared to that 
of g-Gaussian distributed artificial signals. The system size of the 
Ising model is L = 124. We set the same parameters g and q for the 
artificial signals and the model. The data length is 1.6 million. 



the ordered phase, the slope of Fr{n) starts to shift at small 
scales, indicating a multifractal behaviour within these scales. 
Such region is broader when the system is closer to the critical 
point. At the critical point we observe multifractal behaviour 
at all scales we measure. When continuing increasing the tem- 
perature, the system is in the disordered phase. The value of 
the parameter q from the g-Gaussian fit is increasing while 
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FIG. S9: System size dependence of the multifractal behaviour of the 
Ising model. The data length is 1.6 million. 



We then investigate how the multifractal behaviour depends 
on the size of a specific system. Such dependence should be 
manifested from the values of two parameters q and a. To 
see this, we investigate the Ising model with different sizes 
L — 124 and L ~ 512. We first shuffle the data and eliminate 
the dynamics in signals. We then apply the MFDFA and find 
that with the same parameter q, the difference in the observed 
multifractal behaviour is owing to different values of the scale 
parameter a (see Fig. IS9I ). When q is fixed, the larger sys- 
tem with smaller value of a shows multifractal behaviour in 
broader region than the smaller system shows. 



